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Generalized  Gaussian  quadratures  appear  to  have  been  introduced  by  Markov  [11,12]  late 
in  the  last  century,  and  have  been  studied  in  great  detail  as  a  part  of  modern  analysis 
(see  [2,8,9]).  They  have  not  been  widely  used  as  a  computational  tool,  in  part  due  to 
absence  of  effective  numerical  schemes  for  their  construction.  Recently,  a  numerical  scheme 
was  introduced  for  the  design  of  such  quadratures  (see  [10]);  numerical  results  presented 
in  [10]  indicate  that  such  quadratures  dramatically  reduce  the  computational  cost  of  the 
evaluation  of  integrals  under  certain  conditions.  In  this  paper,  we  modify  the  approach  of 
[10],  improving  the  stability  of  the  scheme  and  extending  its  range  of  applicability.  The 
performance  of  the  method  is  illustrated  with  several  numerical  examples. 
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1  Introduction 


Generalized  Gaussian  quadratures  appear  to  have  been  introduced  by  Markov  [11,  12]  late  in 
the  last  century.  More  recent  expositions  include  those  by  Krein  [9]  and  Karlin  [8].  Those 
expositions  contain  proofs  of  the  existence  of  such  quadratures  for  wide  classes  of  functions; 
however  they  do  not  describe  a  numerical  procedure  for  obtaining  the  quadrature  weights  and 
nodes. 

Recently,  a  paper  by  Ma  et  al.  [10]  described  a  numerical  algorithm  for  obtaining  such 
quadratures.  In  [10],  a  version  of  Newton’s  method  is  introduced  for  the  determination  of 
nodes  and  weights  of  generalized  Gaussian  quadratures.  The  procedure  of  [10]  guarantees  the 
convergence  of  the  Newton  algorithm  provided  it  is  started  sufficiently  close  to  the  solution 
(whose  existence  is  proven  in  [11,  9,  8]),  and  utilizes  a  continuation  procedure  to  provide  such 
starting  points.  The  present  paper  describes  a  variation  of  that  algorithm,  which  consists 
mainly  of  two  major  changes.  The  first  change  is  that  an  entirely  different  continuation  scheme 
is  used;  with  the  new  continuation  scheme,  the  algorithm  is  considerably  more  robust.  The 
second  change  is  the  addition  of  a  preprocessing  step  which,  given  as  input  a  large  class  of 
functions,  uses  the  singular  value  decomposition  to  produce  a  set  of  basis  functions  suitable  for 
the  algorithm. 

Since  a  substantial  fraction  of  the  algorithm  is  changed,  this  paper  is  written  as  a  repetition 
of  [10],  rather  than  as  a  list  of  changes;  however,  the  portions  dealing  with  quadratures  for 
functions  with  end-point  singularities  are  omitted. 

This  paper  is  organized  in  the  following  manner.  Section  2  summarizes  the  necessary  ma¬ 
terial  from  [9]  and  [8].  Section  3  briefly  describes  certain  standard  numerical  tools  used  by 
the  algorithm.  Section  4  contains  various  analytical  results  to  be  used  in  the  construction  of 
the  algorithm.  Section  5  describes  the  algorithm  in  detail.  Finally,  Section  6  contains  several 
numerical  examples;  the  actual  nodes  and  weights  obtained  in  Section  6  are  listed  in  Tables 
1-14. 
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2  MathematiccJ  Preliminaries 


We  start  by  introducing  some  notation.  Given  a  finite  sequence  of  real  numbers  xi  <  X2  < 
•  •  •  <  ^ni  let  the  sequence  mi,. .  .,mn  be  defined  as  follows. 


mi 

=  0, 

mj 

=  0 

if  j  >  1  and  Xj  ^  x^-i , 

mj 

=  ;  - 1 

if  j  >  1  and  Xj  =  xj^i  =  . . .  =  xi, 

mj 

=  k 

if  j  >  k  +  1  and  Xj  =  x^^i  =  . . .  =  xj^k  7^  ^j-k-i 

2.1  Chebyshev  systems 


Definition  2.1  A  sequence  of  functions  (f>i,...,4>k  will  be  referred  to  as  a  Chebyshev  system 
on  [a,  b]  if  each  of  them  is  continuous  and  the  determinant 


4>n{xi)  •••  (i)n{Xn) 

is  nonzero  for  any  sequence  of  points  xi,. .  .,Xn  such  that  a  <  Xi  <  X2  . .  ■  <  Xn  <  b. 


(2) 


An  alternate  definition  of  a  Chebyshev  system  is  that  any  linear  combination  of  the  functions 
with  nonzero  coefficients  has  no  more  than  n  zeros. 

A  related  definition  is  that  of  an  extended  Chebyshev  system. 


Definition  2.2  Given  a  set  of  functions  which  are  continuously  differentiable  on  an 

interval  [a,  6],  and  given  a  sequence  of  points  , . . . ,  such  that  a  <  <  X2  <  . . .  <  <  6, 

let  the  sequence  mi,..  .,mn  be  defined  by  (1).  Let  the  matrix  C(xi,...,x„)  =  [cij]  be  defined 
by  the  formula 


c 


ij 


(3) 


in  which  ^^{xj)  is  taken  to  be  the  function  value.  Then  4>i,...,4>n  will  be  referred  to  as  an 
extended  Chebyshev  system  on  [a,  6]  if  the  determinant  |C(a;i, . .  .,Xn)|  is  nonzero  for  all  such 
sequences  Xi . 


Remark  2.1  It  is  obvious  from  Definition  2.2  that  an  extended  Chebyshev  system  is  a  special 
case  of  the  Chebyshev  system.  The  additional  constraint  is  that  the  successive  points  Xi  at  which 
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the  function  is  sampled  to  form  the  matrix  may  be  identical;  in  that  case,  for  each  duplicated 
point,  the  first  corresponding  column  contains  the  function  values,  the  second  column  contains 
the  first  derivatives  of  the  functions,  the  third  column  contains  the  second  derivatives  of  the 
functions,  and  so  forth;  this  matrix  must  also  he  nonsingular. 

Examples  of  Chebyshev  and  extended  Chebyshev  systems  include  the  following  (additional 
examples  can  be  found  in  [8]). 

Example  2.1  The  powers  form  an  extended  Chebyshev  system  on  the  interval 

(-00,  oo). 

Example  2.2  The  exponentials  form  an  extended  Chebyshev  system  for 

any  Ai, . . . ,  >  0  on  the  interval  [0,  oo). 


Example  2.3  The  functions  1,  cosx, sinx,  cos2x,sm2x, . . . ,  cosna:, sinna:  form  a  Chebyshev 
system  on  the  interval  x  6  [0,27r). 

2.2  Generalized  Gaussian  quadratures 


The  quadrature  rules  considered  in  this  paper  are  expressions  of  the  form 

(4) 

j=i 

where  the  points  Xj  €  R  and  coefficients  Wj  €  R  are  referred  to  as  the  nodes  and  weights  of  the 
quadrature,  respectively.  They  serve  as  approximations  to  integrals  of  the  form 


where  a;  is  a  non-negative  function  to  be  referred  to  as  the  weight  function. 

Quadratures  are  typically  chosen  so  that  the  quadrature  (4)  is  equal  to  the  desired  integral 
(5)  for  some  set  of  functions,  commonly  polynomials  of  some  fixed  order.  Of  these,  the  classical 
Gaussian  quadrature  rules  consist  of  n  nodes  and  integrate  polynomdals  of  order  2n  —  1  exactly; 
these  quadratures  are  used  in  this  paper  as  a  numerical  tool  (see  section  3.2).  In  [10],  the  notion 
of  a  Gaussian  quadrature  was  generalized  Jis  follows: 
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Definition  2.Z  A  quadrature  formula  will  be  referred  to  as  Gaussian  with  respect  to  a  set  of 
2n  functions  , <t>2n  •  [o?  6]  ^  R  a,nd  a  weight  function  u  :  [a, 5]  -*  R”^,  if  it  consists  of 

n  weights  and  nodes,  and  integrates  the  functions  <f)i  exactly  with  the  weight  function  u  for  all 
i  =  1,.  ..,2n.  The  weights  and  nodes  of  a  Gaussian  quadrature  will  he  referred  to  as  Gaussian 
weights  and  nodes  respectively. 

The  following  theorem  appears  to  be  due  to  Markov  [11,  12];  proofs  of  it  can  also  be  found 
in  [9]  and  [8]. 

Theorem  2.1  Suppose  that  the  functions  (f>2n  '•  [a, 6]  — *•  R  form  a  Chebyshev  system  on 

[o,  6] .  Suppose  in  addition  that  w  ;  [a,  6]  — >•  R  is  non-negative,  and  is  nonzero  at  more  than  n  - 1 
points  on  [o, 6].  Then  there  exists  a  unique  Gaussian  quadrature  for  <^i, . . . ,<^2n  on  [a, 6]  with 
respect  to  the  weight  function  u.  The  weights  of  this  quadrature  are  positive. 

2.3  Total  positivity 

A  concept  closely  related  to  that  of  an  extended  Chebyshev  system  is  that  of  a  extended  totally 
positive  (ETP)  kernel: 

Definition  2.4  Given  a  function  K  :  [a,  6]  x  [c,  d]  — *•  R  which  is  n  times  continuously  differen¬ 
tiable,  and  given  a  sequence  of  points  xi,...,Xn  such  that  c  <  x\  <  X2  <  . . .  <  x^  <  d,  let  the 
sequence  m\,. ..,  7n„  be  defined  by  (1).  Let  the  functions  4>i,...,4>n  be  defined  by  the  formula 

d'^^K 

=  (6) 

in  which  ^^{xj,t)  is  taken  to  be  the  function  value.  Then  K  will  he  referred  to  as  extended 
totally  positive  if  the  functions  d>i ,  •  •  • ,  d’n  form  an  extended  Chebyshev  system  on  [c,  d]  for  all 
such  sequences  of  Xi. 

Examples  of  ETP  kernels  include  the  following  (additional  examples  can  be  found  in  [8]). 
Example  2.4  The  function  is  extended  totally  positive  for  x,t  £  [0,  oo). 

Example  2.5  The  function  is  extended  totally  positive  for  x,t  £  (-oo,  oo). 
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Example  2.6  The  function  l/(x  +  t)  is  extended  totally  positive  for  x^t  £  (0,  oo) 

A  proof  of  the  following  lemma  can  be  found  in,  for  example,  [8]. 

Lemma  2.2  Suppose  that  K  and  L  are  extended  totally  positive  functions  of  two  variables. 
Then  the  function  M  defined  by  the  formula 

rd 

=  J  K{x,s)L{s^t)ds  (7) 

is  extended  totally  positive.  In  other  words ^  if  the  kernels  of  two  integral  operators  are  extended 
totally  positive^  the  kernel  of  the  product  of  the  two  operators  is  extended  totally  positive. 

The  following  theorem  can  be  found  in  [8,  7]. 

Theorem  2.3  Suppose  that  K  :  [a,  6]  x  [a,  6]  R  is  an  extended  totally  positive  kernel.  Then 
the  first  p  eigenfunctions  of  the  integral  operator  T  :  L^[a,  b]  L^[a,  b]  defined  by  the  formula 

{T(t))(x)=  f  K{x,s)<i){s)ds  (8) 

J  a 

constitute  an  extended  Chebyshev  system^  for  any  p>l. 


3  Numerical  Preliminaries 

3.1  Newton’s  Method 


In  this  section  we  discuss  two  well-known  numerical  techniques:  Newton’s  method  and  the 
continuation  method.  A  more  detailed  discussion  of  these  techniques  can  be  found,  for  example, 
in  [14]. 

Newton’s  method  is  an  iterative  method  for  the  solution  of  nonlinear  systems  of  equations 
of  the  form  F{x)  =  0,  where  F  :  R"  ^  R"  is  a  continuously  differentiable  function  of  the  form 


F{x) 


(  /i(x)  \ 

/2(X) 


\  fn{x)  ) 


(9) 
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and  X  =  (xi,.  ..,Xn)^.  The  method  uses  the  Jacobian  matrix  J  of  F,  which  is  defined  by  the 
formula 

^  •••  §ki^)  ^ 

J{x)  =  :  :  (10) 

V  •••  / 

Lemma  3.1  Suppose  that 

F{y)  =  0  (11) 

with  F  :  ^  R’^  defined  by  (9)^  and  that  \  J(y)\  ^  0,  with  \  J{y)\  denoting  the  determinant  of 

the  matrix  J{x)  defined  in  (10)  at  the  point  y.  Given  a  starting  point  yo  €  R^,  let  the  sequence 
2/1?  2/25  defined  by  the  formula 

yk+i  =  yk-  J~^{yk)F(yk)-  (12) 

Then  there  exists  a  positive  real  numbers  such  that  for  any  yo  satisfying  the  inequality  ||t/o“2/|| 

Sf  the  sequence  (12)  converges  to  y  quadratically,  that  is,  there  exists  a  positive  real  number  a 
such  that 

Wvk+i  -  y\\  <  a\\yk  -  yW^ ■  (13) 

3.1.1  Continuation  method 

In  order  for  Newton’s  method  to  converge,  the  starting  point  which  is  provided  to  it  must  be 
close  to  the  desired  solution.  One  scheme  for  generating  such  starting  points  is  the  continuation 
method,  which  is  as  follows. 

Suppose  that  in  addition  to  the  function  F  :  R^  R^  whose  zero  is  to  be  found,  another 
function  G  :  [0, 1]  x  R’^  ^  R’^  is  available  which  possesses  the  following  properties: 

•  1.  For  any  x  e  R^, 

G(l,x)  =  F(x).  (14) 

•  2.  The  solution  of  the  equation 

G(0,x)=0  (15) 


is  known. 


•  3.  For  all  t  €  [0, 1],  the  equation 


G(t,x)  =  0  (16) 

has  a  unique  solution  x  such  that  the  conditions  of  Lemma  3.1  are  satisfied. 

•  4.  The  solution  x  is  a  continuous  function  of  t. 

If  these  conditions  are  met,  an  algorithm  for  the  solution  of  F{x)  =  0  is  as  follows.  Let  the 
points  ti,ioT  i=  1, . . . ,  m,  be  defined  by  the  formula  f,-  =  i/m.  Solve  in  succession  the  equations 

G{ti,x)  =  0, 

G{t2,x)  =  0, 

G{t^,x)  =  0  (17) 

using  Newton’s  method,  with  the  starting  point  for  Newton’s  method  for  each  equation  taken 
to  be  the  solution  of  the  preceding  equation.  The  solution  x  of  the  final  equation  G{tm,  x)  =  0 
is,  by  (14),  identical  to  the  solution  of  the  desired  equation  F{x)  =  0.  Obviously,  for  sufficiently 
large  m,  Newton’s  method  is  guaranteed  by  Lemma  3.1  to  converge  at  each  step. 

Remark  3.1  In  practice,  it  is  desirable  to  choose  the  smallest  m  for  which  the  above  algorithm 
will  work,  in  order  to  reduce  the  computational  cost  of  the  scheme.  On  the  other  hand,  the 
largest  step  ti  —  for  which  the  Newton  method  will  converge  commonly  varies  as  a  function 
oft.  Thus,  in  this  paper,  we  use  an  adaptive  version  of  the  scheme. 

3.2  Gaussian  integration  and  interpolation 

Classical  Gaussian  quadrature  rules  are  a  well-known  numerical  tool  (see  for  instance  [14]); 
they  integrate  polynomials  of  order  2n  -  1  exactly  with  respect  to  some  weight  function,  and 
consist  of  n  weights  and  nodes.  A  variety  of  Gaussian  quadratures  were  analyzed  in  the  last 
century,  each  being  defined  by  a  distinct  weight  function.  Of  these,  the  algorithm  presented  in 
this  paper  uses  only  the  Gaussian  quadratures  for  the  weight  function  u>{x)  =  1  on  the  region 
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of  integration  [-1, 1].  These  quadratures  are  closely  associated  with  the  Legendre  polynomials; 
we  will  refer  to  their  nodes  as  Legendre  nodes. 

Another  numerical  tool  used  in  this  paper  is  polynomial  interpolation  on  Legendre  nodes. 
Interpolation  refers  to  the  following  problem:  given  two  finite  real  sequences  /i, . . . ,  /n  €  R  and 
x\,...,Xn  €  [c, 6],  construct  a  function  /  :  [a, 6]  ^  R  such  that  /(a:,)  =  /,■  for  all  i  =  1, . . . , n. 
An  interpolation  scheme  is  referred  to  as  linear  if  the  function  /  depends  linearly  on  the  values 
/,.  One  linear  interpolation  scheme  is  polynomial  interpolation,  in  which  the  interpolating 
function  /  is  a  polynomial  of  degree  n  -  1.  As  is  well-known,  such  a  polynomial  always  exists 
and  is  unique.  However,  in  general  two  numerical  difficulties  arise  with  polynomial  interpolation 
using  polynomials  of  high  order.  The  first  is  that  for  many  sequences  of  points  i,-,  the  values 
of  the  interpolating  polynomial  between  the  points  Xi  are  not  well- conditioned  as  a  function 
of  the  values  /,■  to  be  interpolated.  The  second  is  that  even  for  those  sequences  of  points 
where  the  computation  of  the  values  of  the  interpolating  polynomial  is  well-conditioned,  the 
computation  of  the  coefficients  of  the  power  series  of  the  interpolating  polynomial  is  extremely 
ill-conditioned. 

As  is  well-known,  these  difficulties  do  not  arise  if  the  points  x,-  are  taken  to  be  Chebyshev 
nodes  and  the  interpolating  polynomial  is  computed  as  a  series  of  Chebyshev  polynomials  rather 
than  as  a  power  series.  As  the  following  lemma  shows,  the  difficulties  also  do  not  arise  if  the 
points  Xi  are  taken  to  be  Legendre  nodes  and  the  interpolating  polynomial  is  computed  as 
a  series  of  Legendre  polynomials.  The  lemma  makes  use  of  the  following  properties  of  the 
Legendre  polynomials:  first,  that  the  f’th  Legendre  polynomial  Pi  has  degree  i;  second,  that 
the  polynomials  Pi  form  an  orthonormal  system  of  functions  on  [—1,1]. 

Lemma  3.2  Suppose  that  xi,...,x„  €  [—1,1]  the  Legendre  nodes  of  order  n,  and  that 
Wi,...,Wn  €  R  are  the  associated  Gaussian  weights.  Given  a  sequence  /i,...,/n  €  R,  let 
p  :  [— 1,1]  ^  R  6e  the  interpolating  polynomial  of  degree  n  —  1  such  that  p{xi)  =  /,•  for  all 
i  =  1, . . . ,  n,  and  let  cq,  . . . ,  c„_i  be  the  coefficients  of  the  Legendre  series  of  p;  that  is, 

P{x)  =  '^CiPi{x),  (18) 

i=Q 
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(19) 


where  Pi{x)  is  the  iHh  Legendre  polynomial.  Then  the  following  relation  holds: 

Wifi  =  /  Pixfdx  =  £  c?. 

t=l  t=0 

Proof.  The  second  equality  of  (19)  follows  from  (18)  and  the  orthonormality  of  the  Legendre 
polynomials.  The  first  equality  may  be  proven  as  follows:  the  polynomial  p  has  degree  n  —  1, 
thus  its  square  has  degree  2n  —  2.  Since  the  Gaussian  quadrature  integrates  exactly  all  polyno¬ 
mials  up  to  order  2n  —  1,  it  integrates  exactly;  thus  the  first  equality  of  (19)  holds.  □ 

3.3  Singular  value  decomposition 

The  singular  value  decomposition  (SVD)  is  a  ubiquitous  tool  in  numerical  analysis,  which  is 
given  for  the  case  of  real  matrices  by  the  following  lemma  (see,  for  instance,  [3]  for  more  details). 

Lemma  3.3  For  any  n  x  m  real  matrix  A,  there  exists  an  n  X  p  matrix  U  with  orthonormal 
columns,  an  m  x  p  matrix  V  with  orthonormal  columns,  and  a  p  x  p  real  diagonal  matrix 
S  =  [sjj]  whose  diagonal  entries  are  non-negative,  such  that  A  =  USV*  and  that  sa  <  Sj+i,,+i 
for  all  i  =  1, . .  .,p  —  1. 

The  diagonal  entries  s,-  of  S  are  called  singular  values;  the  columns  of  the  matrix  V  are 
called  right  singular  vectors;  the  columns  of  the  matrix  U  are  called  left  singular  vectors. 

3.4  Singular  value  decomposition  of  integral  operators 

This  section,  which  follows  [5],  contains  an  existence  theorem  for  a  factorization  of  integral 
operators.  The  operators  T  :  L^[c,  d]  —>■  L^[a,  6]  to  which  it  applies  are  of  the  form 

{Tf){x)  =  j\{x,t)f{t)dt.  (20) 

in  which  the  function  K  :  [a,  6]  x  [c,(i]  —»•  R  is  referred  to  as  the  kernel  of  the  operator  T. 
Throughout  this  section,  it  will  be  assumed  that  all  functions  are  square-integrable;  the  term 
“norm”  will  mean  the  norm. 

The  following  theorem,  which  defines  the  factorization,  is  proven  in  a  more  general  form  as 
Theorem  VI.  17  in  [13]. 
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Theorem  3.4  Suppose  that  the  function  K  :  [a,  6]  x[c^d\ 'R.  is  square  integrable.  Then  there 
exist  two  orthonormal  sequences  of  functions  Ui  :  [a,  6]  — >  R  and  Vi  :  [c,  d]  R  and  a  sequence 
Si  €  R,  /or  i  =  1, . . .  ,  cx),  such  that 

CO 

K{x,t)  =  '^Ui{x)siVi{t)  (21) 

t=l 

and  that  >  52  >  . . .  >  0.  The  sequence  Si  is  uniquely  determined  by  K,  Furthermore,  the 
functions  Vi  are  eigenfunctions  of  the  operator  T*T,  where  T  is  defined  by  (20),  and  the  values 
Si  are  the  square  roots  of  the  eigenvalues  ofT*T. 


By  analogy  to  the  finite- dimensional  case,  we  will  refer  to  this  factorization  as  the  singular 
value  decomposition.  We  will  refer  to  the  functions  Ui  as  left  singular  functions  of  K  (or  of  T), 
to  Vi  as  right  singular  functions,  and  to  Si  as  singular  values. 

As  is  the  case  for  the  discrete  singular  value  decomposition,  this  decomposition  can  be  used 
to  construct  an  approximation  to  the  function  by  discarding  small  singular  values  and  the 
associated  singular  functions: 

K{x,t)  ~'^Ui{x)$iVi{t).  (22) 

t=i 

The  error  of  this  approximation  can  then  be  computed  from  (21); 

P  OO 

K{x,t)-'^Ui{x)siVi{t)=  Ui{x)siVi{t),  (23) 

Z=1  2=:p-|-l 


and,  therefore. 


K{x,  0  -  H  tii(x)5,Vi(t) 


2=1 


2=p+l 


Using  (24),  we  will  be  approximating  integrals 

f  K{x,t)u(x)dx 
Ja 


(24) 


(25) 


by  the  formula 


r  r  ^ 

/  K{x,t)u{x)dx  ^  Ui{x)siVi{t)u{x)dx 

Ja  Ja 

V  .h 

-  /  Ui{x)(jj{x)dx, 

2  =  1 


(26) 
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Thus,  a  quadrature  which  is  exact  for  each  of  the  integrals 


fb 

/  Ui{x)cj(x)dx^  (27) 

Ja 

for  i  =  1, . .  .,,p,  is  an  approximate  quadrature  for  integrals  of  the  form  (25). 

The  following  theorem  shows  that  if  an  operator  is  extended  totally  positive,  its  singular 
functions  form  an  extended  Chebyshev  system. 


Theorem  3.5  Suppose  that  K  :  [a,  6]  x  [c,  d]  — »•  R  is  extended  totally  positive.  Then  the  first  p 
left  singular  functions  of  K  form  an  extended  Chebyshev  system,  for  any  p;  likewise  the  first  p 
right  singular  functions  of  K  form  an  extended  Chebyshev  system,  for  any  p. 


Proof.  Let  the  integral  operator  T  :  L^[c,d]  L^[a,b]  be  defined  by  the  formula 

{Tf){x)  =  j\{x,t)f{t)dt,  (28) 

and  the  function  L  :  [a,  b]  — »■  [a,  6]  be  defined  by  the  formula 

L{x,t)  =  J  K{x,s)K{t,s)ds.  (29) 

Clearly,  the  integral  operator  5  :  L^[o,  6]  L‘^[a,b]  defined  by  the  formula  5  =  T*T  has  the 

kernel  L: 


{S4>){x) 


nK  (x,  s)K  (t,  s)ds(f>(t)dt 

_ 

I  L{x,t)4>{t)dt. 

Ja 


(30) 


Since  K  is  extended  totally  positive,  due  to  Lemma  2.2,  L  is  also  extended  totally  positive. 
Thus,  by  Theorem  2.3,  the  eigenfunctions  of  S  constitute  an  extended  Chebyshev  system.  By 
Theorem  3.4,  these  eigenfunctions  are  identical  to  the  left  singular  functions  of  T,  which  proves 
that  the  first  p  left  singular  functions  of  T  constitute  an  extended  Chebyshev  system,  for  any 
p.  The  proof  for  the  right  singular  functions  is  identical.  □ 
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4  Analytical  Apparatus 

4.1  Convergence  of  Newton’s  method 


In  this  section,  we  observe  that  the  nodes  and  the  weights  of  a  Gaussian  quadrature  satisfy  a 
certain  system  of  nonlinear  equations.  We  then  prove  that  the  Newton  method  for  this  system  of 
equations  is  always  quadraticaJly  convergent,  provided  the  functions  to  be  integrated  constitute 
an  extended  Chebyshev  system. 

Given  a  set  of  functions  and  a  weight  function  oj,  the  Gaussian  quadrature  is 

defined  by  the  system  of  equations 


j=i 

n 

^  Wi4>2{Xi) 


i=i 


f  <j>i{x)u{x)dx, 
Ja 

fb 

/  4>2ix)u}{x)dx, 
Ja 


f  4>2n{x)ojix)da 

Ja 


(31) 


(see  Definition  2.3).  Let  the  left  hand  sides  of  these  equations  be  denoted  by  fi  through  /2n- 
Then  each  /,  is  a  function  of  the  weights  wi,. .  .,Wn  and  nodes  xi , . . . , of  the  quadrature. 
Its  partial  derivatives  are  given  by  the  obvious  formulae 


dh 

dwi 


4>k{Xi)., 


dfk 

dxi 


Wi(i>'k{xi). 


Thus,  the  Jacobian  matrix  of  the  system  (31)  is 


(32) 

(33) 


J(Xi,.  .  .,Xn,U)i, 


/  <^i(a:i) 


Wn)  = 


4>\{Xn)  Wi(?i>i(xi) 


.  (34) 


V  ^27i(^i)  '**  ^2n(^n)  **’  "^71^271(^^1)  / 

Lemma  4*1  Suppose  that  the  functions  . .  .,<^2n  form  an  extended  Chebyshev  system.  Let 
the  Gaussian  quadrature  for  these  functions  be  denoted  by  Wi  and  x^.  Then  the  determi¬ 
nant  of  J  is  nonzero  at  the  point  which  constitutes  the  Gaussian  quadrature;  in  other  words y 
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Proof.  It  is  immediately  obvious  from  (34)  that 


Wi  •  'W2 


Wn-\  •  Wn 


(35) 


*^271(2^1)  ■  ■  ■  4*2n{^n)  ^2n(®l)  ‘  ’  '  ^2n(^n) 

If  ,  <j)2n  form  an  extended  Chebyshev  system,  then  by  Theorem  2.1 ,  the  weights  , . . . ,  ti)„ 
of  the  Gaussian  quadrature  are  positive.  In  addition,  by  the  definition  of  an  extended  Chebyshev 
system,  the  determinant  in  the  right  hand  side  of  (35)  is  nonzero.  Thus 


\J{xu...,Xn,Wi,...,Wi)\  7^  0. 


(36) 


□ 


Corollary  4.2  Under  the  conditions  of  Lemma  4-U  ^he  Gaussian  weights  and  nodes  depend 
continuously  on  the  weight  function. 

4.2  Linear  interpolation 

Given  a  collection  of  n  points  xi,X2, .  ■  ■,Xn  €  an  interpolation  scheme  with  the  nodes 

xi,X2,  ■ .  .,Xn  is  a  linear  mapping  T  :  R”  — >•  L°°[a,  6]  such  that  for  any  point  j/  6  R”, 

Liy){xi)  =  yi,  (37) 

for  all  i  =  1,2,  For  a  function  /  €  6],  the  error  5l(/)  of  interpolation  of  the 

interpolation  scheme  L  is  defined  by  the  formula 

h{f)  =  max  \  fix)  -  T(T’)(a:)|  (38) 

xG[a,o] 

where  F  =  (/(xi), /(X2), . . /(xn))^- 

The  following  lemma  serves  as  a  justification  for  the  notation  we  use  in  Section  4.3  for  linear 
interpolation  schemes. 


13 


Lemma  4.3  Suppose  L  :  R”  — L^[a^b]  is  a  linear  interpolation  scheme  with  n  nodes 
xi,...,Xn  €  [a,  6].  Then  there  exists  a  sequence  of  functions  ai,...,o:n  :  [c,  6]  -*  R  such 
that  for  any  vector  f  6  R",  with  elements  /  =  (/i, . . . ,  fn)^, 

{Lf){x)  =  f2fMx),  (39) 

t=i 

for  all  X  €  [a,  ^>]- 

Proof.  Let  the  vectors  Ci, . . . ,  Cn  €  R"  with  elements  e,  =  (e,i, . . . ,  e,n)^  be  the  standard  basis 
in  R";  that  is,  e,,-  =  1  for  all  i  =  1, . . . ,  n,  and  ej^  =  0  for  aU  t,  j  =  1, . . . ,  n  such  that  i  ^  j.  Let 
the  functions  ai, . . . ,  :  [a,  6]  — >■  R  be  defined  by  the  formula  ai  =  Xe,-.  Since  the  interpolation 

scheme  L  is  linear,  for  any  vector  /  €  R”  with  elements  /  =  {fi,  •  •  •■,  fn)^ ,  for  any  point 

X  €  [a,  &], 

(L/)(i)  =  (x(|:/iei))  (x) 

=  l^/i(M)(x) 

t=l 

□ 

In  the  case  of  polynomial  interpolation,  the  functions  ai  are  referred  to  as  Lagrange  polynomials; 
by  analogy  to  that  case,  we  will  in  general  refer  to  the  functions  as  the  Lagrange  functions 
of  the  interpolation  scheme. 

The  following  lemma  provides  an  error  bound  for  approximation  of  a  function  of  two  vari¬ 
ables  using  two  one- dimensional  interpolation  formulae,  expressed  in  terms  of  error  bounds 
for  each  one-dimensional  interpolation  scheme  applied  separately.  Its  proof  is  an  exercise  in 
elementary  analysis,  and  is  omitted. 

Lemma  4.4  Suppose  that  a:i,X2, . .  .,a:n  €  [a, 6]  and  €  [c,d]  are  two  finite  real 

sequences,  and  that  Cki,a2,  •  •  *,otn  •  [a, 6]  — >  R  and  /?i, /?29  •  •  •  ? ^  [c? ^  ore  two  sequences 
of  bounded  functions.  Suppose  further  that  Lj  :  R^  ^  L^[a,  6]  is  an  interpolation  formula 
with  the  nodes  xi,...,X7i  and  Lagrange  functions  ai,...,a7i,  and  L2  :  R^  ^  L°®[c, d]  is  an 
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interpolation  formula  with  the  nodes  ct'^^d  Lagrange  functions  Suppose 

that  Tf  is  such  that 

^lai(x)|<7/,  (41) 

t=i 

for  all  X  €  [a,  6],  and 

m 

i=i 

for  all  t  £  [c,d].  Finally,  suppose  that  K  is  a  function  [a,  6]  x  [c,d]  -*  R,  and  that  for  all 
X  £  [a,  6]  and  t  £  [c,  d], 


K{x,t)-J2K{xi,t)ai{x) 


i=l 


and 


Then 


K{x,t)-J2K{x,t^)0,{t) 

i=i 


K{x,t)-J2Y1 

t=l  j  =  l 


<  £ 


<  £. 


<  <1  +-^), 


(43) 

(44) 

(45) 


for  all  X  £  [a,  6]  and  t  £  [c,  d]. 

4.3  Approximation  of  SVD  of  an  integral  operator 

This  section  describes  a  numerical  procedure  for  computing  an  approximation  to  the  singular 
value  decomposition  of  an  integral  operator. 

The  algorithm  uses  quadratures  which  possess  the  following  property. 

Definition  4.1  We  will  say  that  the  combination  of  a  quadrature  and  a  linear  interpolation 
scheme  preserves  inner  products  on  an  interval  [o,  6]  if  it  possesses  the  following  properties. 

•  1.  The  nodes  of  the  quadrature  are  identical  to  the  nodes  of  the  interpolation  scheme. 

•  2.  The  quadrature  integrates  exactly  any  product  of  two  interpolated  functions;  that  is, 
for  any  two  functions  /,  5  :  [a,  6]  ^  R  produced  by  the  interpolation  scheme,  the  integral 

rb 


/ 


f{x)g{x)dx 


(46) 


is  computed  exactly  by  the  quadrature. 
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Quadratures  and  interpolation  schemes  which  possess  this  property  include: 

Example  4.1  The  combination  of  a  (classical)  Gaussian  quadrature  at  Legendre  nodes  and 
polynomial  interpolation  at  the  same  nodes  preserves  inner  products^  since  polynomial  interpo¬ 
lation  on  n  nodes  produces  an  interpolating  polynomial  of  order  n  —  1,  the  product  of  two  such 
polynomials  is  a  polynomial  of  order  2n  —  2,  and  a  Gaussian  quadrature  integrates  exactly  all 
polynomials  up  to  order  2n  -  I, 

Example  4.2  If  an  interval  is  broken  into  several  subintervals,  and  a  quadrature  and  interpo¬ 
lation  scheme  which  preserves  inner  products  is  used  on  each  subintervalj  then  the  arrangement 
as  a  whole  preserves  inner  products  on  the  original  interval.  (This  follows  directly  from  the 
definition.) 

Example  4.3  The  combination  of  the  trapezoidal  rule  on  the  interval  [0,27r],  and  Fourier  in¬ 
terpolation  (using  the  interpolation  functions  1,  cos  x,  sin  x,  cos  2x,  sin  2a:, ,  cos  nx,  sin  nx )  pre¬ 
serves  inner  products. 

The  algorithm  takes  as  input  a  function  K  :  [a,  b]  x  [c,  ci]  — >  R.  It  uses  the  following  numerical 
tools: 

•  1.  A  quadrature  and  a  linear  interpolation  scheme  on  the  interval  [a,  6]  which  preserve 
inner  products.  Let  the  weights  and  nodes  of  this  quadrature  be  denoted  by  wf., . . . ,  G 
R  and  xi, . . €  [a,  6]  respectively.  Let  the  Lagrange  functions  (see  Section  4.2)  of  the 
interpolation  scheme  be  denoted  by  ai an  :  [a,  6]  ^  R. 

•  2.  A  quadrature  and  a  linear  interpolation  scheme  on  the  interval  [c,d]  which  preserve 

inner  products.  Let  the  weights  and  nodes  of  this  quadrature  be  denoted  by  . . . ,  € 

R  and  €  [c,d]  respectively.  Let  the  Lagrange  functions  of  the  interpolation 

scheme  be  denoted  by  /3i, . . . ,  :  [c?  d]  R- 

As  will  be  shown  below,  the  accuracy  of  the  algorithm  is  then  determined  by  the  accuracy  to 
which  the  above  two  interpolation  schemes  approximate  K. 
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The  output  of  the  algorithm  is  a  sequence  of  functions  tii , . . . ,  Up  :  [a,  6]  R,  a  sequence  of 
functions  ui, . .  .,t;p  :  [c,  d]  ^  R,  and  a  sequence  of  singular  values  5i,...,Sp  £  R,  which  form 
an  approximation  to  the  singular  value  decomposition  of  K. 

Description  of  the  algorithm: 

•  1.  Construct  the  n  x  m  matrix  A  =  [oy]  defined  by  the  formula 

c,j  =  K{xi,  tj )  ^wf-w).  (47) 

•  2.  Compute  the  singular  value  decomposition  of  A,  to  produce  the  factorization 

A  =  USV\  (48) 

where  U  =  [ttij]  is  an  n  x  p  matrix  with  orthonormal  columns,  V  =  [vij]  is  an  m  x  p 
matrix  with  orthonormal  columns,  and  5  is  a  p  X  p  diagonal  matrix  whose  j’th  diagonal 
entry  is  Sj . 

•  3.  Construct  the  n  x  p  matrix  U  =  [uij]  and  the  m  x  p  matrix  V  =  [vij]  defined  by  the 
formulae 


Uik  =  Uik/yJ^, 

Vjk  =  Vjkl\Jw]. 

(49) 

(50) 

•  4.  For  any  points  x  €  [o,6]  and  t  €  [c,d],  evaluate  the  functions 

Ufc  :  [c,  d]  — »■  R  via  the  formulae 

Uk  :  [a,b\  — >•  R  and 

n 

‘^k(x)  =  '^Uik  •  Oii{x), 
i=l 

(51) 

m 

■  I3jit), 
i=i 

(52) 

for  all  A:  =  1, . .  .,p. 

Theorem  4.5  Suppose  that  the  combination  of  the  quadrature  with  weights  and  nodes 

, . . . ,  G  R  and  Xi, . . . ,  Xn  €  [a,b],  respectively,  and  the  interpolation  scheme  with  Lagrange 

functions  , . . . ,  :  [a,  6]  — »  R,  preserves  inner  products  on  [a,  6]. 
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Suppose  in  addition  that  the  combination  of  the  quadrature  with  weights  and  nodes 
Wj, . . . ,  €  R  and  ti, . . . ,  €  [c,  ci],  respectively,  and  the  interpolation  scheme  with  Lagrange 

functions  ^i, . . . ,  Sm  '•  [c,  d]  R,  preserves  inner  products  on  [c,  d]. 

For  any  function  K  :  [o,  6]  x  [c,  d]  — >■  R,  let  u,-  :  [o,  6]  K,  Vi  :  [c,  d]  R,  and  Si  £  Ji  be 
defined  in  (47)-(52),  for  all  i  =  1, . .  .,p.  Then 

•  1.  The  functions  Ui  are  orthonormal,  i.e. 

f  Ui{x)uk{x)dx  =  Sik  (53) 

Ja 

for  all  z,  =  1, . .  with  Sik  I'he  Kronecker  symbol  (Sij  —  \  if  i  0  otherwise). 


•  2.  The  functions  Vi  are  orthonormal,  i.e 

rd 


)  Vi{t)vk{t)dx  =  Sik 


(54) 


for  all  i,k  =  1,. .  .,p. 


3.  The  function  K  :  [a,  b]  x  [c,  (i]  R  defined  by  the  formula 

p 

K{x,t)  =  J2sjUj(x)vj{t), 
j=l 


(55) 


is  identical  to  the  function  produced  by  sampling  K  on  the  grid  of  points  then 

interpolating  with  the  two  interpolation  schemes.  That  is, 


jk{x,t)  =  ^^ii:(a:i,ij)Q!i(a:)/3j(i). 
t=i  i=i 

Proof.  We  first  prove  (56).  Combining  (51),  (52),  and  (55),  we  have 

P  /  n  \  ^  ^ 

K{x,t)  =  ^Ufc(x,>.(2)) 

k=i  \»=i  /  \i=i  ) 

n  m  /  p  \ 

=  EE  Xi)skVkiw^}  ]  ai{x)/3j{t) 

t=l  j=l  \k^l  / 

n  m  /  p  _  _  \ 

t=ij=i  \k=i  J 


(56) 


18 


=  E  (  S  ^ikSkVjkfy/wfw)  ]  ai(x)pj{t) 

i=lj=zl  \A:=1  / 

n  m 

=  £  (°«i/ Q:.(*)/J;(<)- 

t=l  j=l 


(57) 


Now,  (56)  follows  from  the  combination  of  (57)  and  (47). 

We  now  demonstrate  the  orthonormaJity  of  the  functions  u,-.  Since  these  are  functions 
produced  by  interpolation,  and  since  the  quadrature  on  [o,  5]  is  assumed  to  integrate  exactly 
aU  products  of  pairs  of  interpolated  functions, 

fb 

/  Ui{x)ukix)dx  =  '^WjUi{xj)uk{xj) 
i=i 
n 


i=i 

n 

i=i 


(58) 


Since  the  last  sum  in  (58)  is  the  inner  product  of  two  columns  of  the  orthonormal  matrix  U 
(see  (48), 

/  ui{x)uk{x)dx  =  Sik-  (59) 

J  a 

The  orthonormality  of  the  functions  Vi  is  proven  in  the  same  manner.  □ 


Remark  4.1  Obviously y  the  above  proof  approximates  the  singular  value  decomposition  of  the 
operator  T  :  i^[c,  d]  — >  Z/^[a,6]  with  the  kernel  K  by  constructing  an  approximation  T  with 
kernel  K  to  the  operator  T  that  is  of  finite  rank,  and  constructing  the  exact  singular  value 
decomposition  of  the  latter. 

Observation  4.2  In  the  preceding  proof,  the  assumption  that  each  combination  of  quadrature 
and  interpolation  scheme  preserves  inner  products  was  used  only  to  demonstrate  the  orthonor- 
mality  of  the  corresponding  singular  functions.  Thus,  if  the  conditions  of  Theorem  4-5  hold, 
with  the  exception  that  the  quadrature  on  [a,b]  does  not  preserve  inner  products,  then  (54)  o>Tid 
(56)  hold  (but,  in  general,  (53)  does  not). 
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Remark  4*3  Theorem  and  Lemma  4-4  generalize  trivially  to  higher  dimensions.  One- 
dimensional  quadratures  and  interpolation  formulae  have  to  be  replaced  with  their  multidimen¬ 
sional  counterparts;  otherwise,  the  proofs  are  unchanged. 

5  Numericad  Algorithm 

This  section  describes  a  numerical  algorithm  for  the  evaluation  of  nodes  and  weights  of  gen¬ 
eralized  Gaussian  quadratures.  The  algorithm’s  input  is  a  sequence  of  functions  • 

[a,  6]  — >•  R  which  form  an  extended  Chebyshev  system  on  [o,  6],  and  a  weight  function  : 
[a,b]  R"*".  Its  output  is  the  weights  and  nodes  of  the  quadrature.  The  main  components  of 

the  algorithm  are  as  follows  (not  listed  in  order  of  execution): 

•  1.  Newton’s  method  is  used  to  solve  (31)  which  defines  the  Gaussian  quadrature. 

•  2.  An  adaptive  version  of  the  continuation  method  (Section  3.1.1)  is  used  to  provide 
starting  points  for  Newton’s  method.  The  continuation  scheme  used  here  is  different  from 
that  used  in  [10];  the  details  of  the  continuation  scheme  and  of  the  method  of  adaption 
are  described  below. 

•  3.  The  algorithm  of  section  4.3  can  be  used  as  an  optional  preprocessing  step,  which  takes 
as  input  a  kernel  of  an  integral  operator  and  produces  its  singular  functions.  The  first  2ra 
of  the  left  singular  functions  are  then  used  as  input  to  the  main  algorithm. 

5.1  Continuation  Scheme 

The  continuation  scheme  used  is  as  follows.  Let  the  weight  functions  u  :  [0, 1]  x  [a,  6]  — »•  R"*"  be 
defined  by  the  formula 

n 

u;(a,x)  =  awi(x)  +  (1  -  a)^6(x  -  Cj),  (60) 

i=i 

where  is  the  weight  function  for  which  a  Gaussian  quadrature  is  desired,  6  denotes  the  Dirac 
delta  function,  and  the  points  Cj  £  [a,  6]  are  arbitrary  distinct  points.  These  weight  functions 
have  the  following  properties: 
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•  1.  With  a  =  1,  the  weight  function  is  equal  to  the  desired  weight  function  a;i,  due  to 
(60). 


•  2.  With  a  =  0,  the  Gaussian  weights  and  nodes  are 

Wj  =  1,  (61) 

Xj  =  Cy,  (02) 

for  j  =  1, . .  .,n,  whatever  the  functions  4>i  are  (since  a;(0,a:)  =  0,  unless  x  =  Cj  for  some 
i  e  [l,n]). 

•  3.  The  quadrature  weights  and  nodes  depend  continuously  on  a  (by  CoroUary  4.2). 

The  intermediate  problems  which  the  continuation  method  solves  are  the  Gaussian  quadratures 
relative  to  the  weight  functions  a;(a,*).  The  scheme  starts  by  setting  a  =  0,  then  increases  a 
in  an  adaptive  manner  until  a  =  1,  as  follows.  A  current  step  size  is  maintained,  by  which  a 
is  incremented  after  each  successful  termination  of  Newton’s  method.  After  each  unsuccessful 
termination  of  Newton’s  method,  the  step  size  is  halved  and  the  algorithm  restarts  from  the 
point  yielded  by  the  last  successful  termination.  After  a  certain  number  of  successful  steps,  the 
current  step  size  is  doubled.  (Experimentally,  the  current  problem  was  found  to  be  well  suited 
to  an  aggressive  mode  of  adaption:  in  the  authors’  implementation,  the  initial  value  of  the  step 
size  was  chosen  to  be  one,  and  the  step  size  was  doubled  after  a  single  successful  termination 
of  Newton’s  method.) 

5.1.1  Starting  points 

The  choice  of  the  points  Cj  was  left  indefinite  above.  In  exact  arithmetic  the  algorithm  would 
converge  for  any  choice  of  distinct  points  (see  Lemma  4.1).  However  the  number  of  steps  of  the 
continuation  method,  and  thus  the  speed  of  execution,  is  affected  by  the  choice.  More  impor¬ 
tantly,  the  numerical  stability  of  the  scheme  might  be  compromised  due  to  poor  conditioning  of 
the  matrix  J  (see  (34)).  Indeed,  while  Lemma  4.1  guarantees  that  the  matrix  J  is  non- singular, 
it  says  nothing  about  its  condition  number.  Thus,  in  the  authors’  implementation,  the  points 
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Cj  used  for  the  production  of  the  quadrature  of  order  n  were  computed  from  the  nodes  Xj  of 
the  quadrature  of  order  n  -  1,  by  the  formulae 


Cl  =:  Xi, 

(63) 

Ci  =  1  H“  ^t)/2, 

II 

.,n-  1, 

(64) 

Cji  = 

(65) 

With  this  choice,  no  failures  to  converge  have  been  encountered  in  the  authors’  experience. 

6  Numerical  examples 

A  variety  of  quadratures  were  generated  to  illustrate  the  performance  of  the  above  algorithm. 
In  each  case  the  preprocessing  step  of  producing  singular  functions  was  used.  This  step  requires 
two  sets  of  quadratures  and  interpolation  schemes,  which  must  approximate  the  desired  kernel 
to  the  desired  accuracy.  These  quadratures  and  interpolation  schemes  were  chosen  so  that 
the  approximation  was  accurate  to  about  the  precision  of  the  arithmetic  that  was  used.  The 
following  combination  of  quadrature  and  interpolation  scheme  which  preserves  inner  products 
was  used:  the  interval  of  integration  was  divided  into  several  subintervals,  and  a  combination 
of  a  (classical)  Gaussian  quadrature  at  Legendre  nodes  and  polynomial  interpolation  was  used 
on  each  subinterval. 

In  each  of  the  following  examples,  the  calculations  were  done  in  extended  precision  (Fortran 
REAL*  16)  arithmetic,  with  the  exception  of  the  last  example,  which  was  done  in  double  precision 
(REAL*8)  arithmetic. 

6.1  Exponentials 

In  this  example  we  construct  quadratures  for  the  integral 

TOO 

/  (66) 

Jo 

under  the  condition  that  1  <  t  <  500.  In  this  case,  the  corresponding  kernel  K  :  [0,oo)  x 
[1,500]  -»■  R  is  given  by 

K{x,t)  =  e-^\  (67) 
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and  is  extended  totally  positive;  thus  its  singular  functions  form  an  extended  Chebyshev  system. 
A  sample  of  the  quadratures  produced  by  the  algorithm  is  included  in  Tables  1-3;  for  double 
precision  accuracy,  a  27-poiEt  quadrature  is  required. 

6.2  Complex  Exponentials 

Here,  we  design  quadratures  for  a  new  version  [5]  of  the  two-dimensional  Fast  Multipole  Method. 
These  quadratures  are  for  the  integral 

r  e-^^dx,  (68) 

Jo 

under  the  condition  that  z  €  C  is  constrained  to  lie  in  the  region  D  of  the  complex  plane  which 
consists  of  the  rectangle  [1,4]  X  [-4,4]  with  a  1  X  1  square  deleted  from  each  of  its  two  left  hand 
corners,  as  depicted  in  Figure  1.  Since  both  the  true  integral  (equal  to  1/z)  and  the  quadrature 
which  approximates  the  integral  are  complex  analytic  on  that  region,  due  to  the  maximum 
modulus  principle  the  maximum  error  of  the  quadrature  is  achieved  on  the  boundary  SD  oi  the 
region.  Accordingly,  the  kernel  whose  singular  functions  were  computed  was  K{x,z)  = 
with  z  varying  over  SD.  A  brief  examination  of  the  resulting  singular  functions  shows  that  they 
do  not  form  a  Chebyshev  system;  if  they  did  so,  the  i’th  function  would  have  i  —  1  zeros,  yet  it 
has  many  more.  Thus  the  algorithm  is  not  guaranteed  to  work;  however,  it  did  so.  A  sample 
of  the  resulting  quadratures  is  included  in  Tables  4-6;  in  this  case  a  quadrature  yielding  double 
precision  accuracy  contains  32  nodes. 

6.3  Exponentials  multiplied  by  Iq 

In  this  example,  quadrature  formulae  are  constructed  for  integrals  of  the  form 

AOO 

/  h{xy)e-^^dx,  (69) 

Jo 

under  the  condition  that  t  €  [1, 500]  and  j/  €  [0,t  —  1];  these  formulae  were  designed  to  be  used 
in  a  version  of  the  one-dimensional  Fast  Multipole  Method  which  is  used  in  an  algorithm  [6]  for 
the  fast  Hankel  transform.  In  this  case  the  singular  functions  produced  by  the  precomputation 
stage  were  extremely  similar  to  those  for  exponentials  alone;  unlike  in  the  case  of  complex 
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Figure  1:  Coefficient  z  of  complex  exponentials  to  be  integrated 

exponentials,  it  is  possible  that  they  form  a  Chebyshev  system.  In  any  case,  the  algorithm 
converged,  producing  a  quadrature  which  required  two  more  nodes  for  double  precision  accuracy 
than  were  required  for  the  integration  of  exponentials  alone  (i.e.  29  nodes).  A  sample  of  the 
resulting  quadratures  is  included  in  Tables  7-9. 

6.4  Exponentials  multiplied  by  Jo 
Here,  we  construct  quadratures  for  the  integral 

roo 

/  Jo{xy)€^'^^dx,  (70) 

under  the  conditions  that  t  G  [1,4]  and  y  €  [0, 4v^],  where  Jo  denotes  the  Bessel  function  of  the 
first  kind  of  order  zero.  These  quadratures  are  used  in  a  new  version  [4]  of  the  three-dimensional 
Fast  Multipole  Method.  Jo  is  given  by  the  well-known  (see  for  instance  [1])  formula 

Jo{y)  =  “  /  cos{y  cos  9)d6,  (71) 

TT  Jo 
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Substituting  (71)  into  (70)  yields  tbe  integral 


cos{xy  cos  6)d6^  e 

I  rir  foo 

=  —  y  J  cos{xy  COS  9)e~^dxd0. 
Thus  a  quadrature  accurate  for  the  integral 

TOO 

/  cos{xy)e~^^dx^ 

Jo 


(72) 


(73) 


under  the  conditions  that  t  €  [1,4]  and  y  G  [0,4v^],  is  also  accurate  for  the  integral  (70)  under 
the  same  conditions  on  y  and  t.  Since  the  function  cos(a:y)e~^*  is  a  harmonic  function  of  y 
and  t,  by  the  maximum  modulus  principle  the  maximum  error  of  a  quadrature  for  (73)  lies 
on  the  boundary  6D  of  the  rectangular  region  t  G  [1,4],  y  G  [0,4\/2].  Accordingly,  the  kernel 
whose  singular  functions  were  computed  was  K{x^z)  =  cos(a;y)e“'^*,  with  (t^y)  varying  over 
SD,  As  in  the  case  of  complex  exponentials,  the  singular  functions  have  too  many  zeros  to  form 
a  Chebyshev  system,  however  the  algorithm  converged.  A  sample  of  the  resulting  quadratures 
is  included  in  Tables  10”14;  for  single  precision  accuracy  22  nodes  are  required. 


6.5  Numerical  Observations 


The  following  observations  were  made  in  the  course  of  our  numerical  experiments. 


•  1.  The  number  of  continuation  steps  required  is  highly  variable;  in  many  cases,  only  one 
step  suiRced  to  produce  the  quadrature;  less  frequently,  up  to  fifty  or  so  continuation 
steps  were  required.  This  variability  occurred  even  between  quadratures  for  successive 
numbers  n  of  nodes,  with  the  same  weight  function  and  kernel  K, 


•  2.  The  algorithm  worked  in  the  cases  where  Theorem  2.1  applied,  and  also  in  cases  where 
it  did  not.  In  the  latter  cases,  it  is  conceivable  that  the  resulting  quadratures  would  have 
negative  weights,  or  that  they  would  not  be  unique.  However,  all  computed  weights  were 
positive,  and,  while  no  systematic  attempt  was  made  to  look  for  non-uniqueness  of  the 
quadratures,  no  instance  of  it  was  observed. 
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Generalizations  and  Applications 


•  1.  The  success  of  the  algorithm  in  instances  where  Theorem  2.1  does  not  apply  sug¬ 
gests  that  further  theoretical  investigation  of  conditions  for  the  existence  of  generalized 
Gaussian  quadratures  would  be  profitable. 

•  2.  An  obvious  generalization  of  these  results  is  to  quadratures  for  integrals  in  more  than 
one  dimension.  However,  such  an  extension  does  not  seem  to  have  been  explored  classi¬ 
cally;  the  authors  are  investigating  a  generalization  of  Theorem  2.1  for  multidimensional 
quadratures. 

•  3.  An  obvious  application  of  the  algorithm  of  this  paper  is  for  the  efficient  evaluation  of 
functions  represented  by  their  integral  transforms  (see  Sections  6.1,  6.2,  6.3,  6.4  above, 
as  well  as  [5]  and  [4]).  The  method  of  steepest  descent  in  the  numerical  complex  analysis 
provides  a  wide  field  of  applications  for  such  algorithms.  . 

•  4.  An  entirely  different  field  of  applications  involves  the  numerical  solution  of  integral 
equations  with  singular  kernels;  of  particular  interest  are  boundary  integral  equations 
of  scattering  theory  on  regions  with  corners.  The  authors  are  currently  pursuing  this 
direction  of  research. 
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Table  1:  Quadratures  for  exponentials 


Quadratures  for  the  integral 

Jroo 

'  €~^*dx, 

0 

under  the  condition  that  1  <  t  <  500. 


N 

Nodes  (x,) 

Weights  (wi) 

Error 

“T 

0.2934661296034111E-02 

0.2122706574797170E-01 

0.8809516681098265E-01 

0.3048205241689060E+00 

0.9407821348001514E+00 

0.2710823671 107057E+01 

0.8078894059616301E-02 

0.3337852721645502E-01 

0.1157432569817795E+00 

0.3589923073929015E+00 

0.1018252445219498E+01 

0.2863049428178813E+01 

0.827E-03 

8 

0.2027451178542047E-02 

0.1244627909236754E-01 

0.4102057941602644E-01 

0 . 1 14493  7922230447E+00 
0.2954175879426304E+00 
0.7210624530246545E+00 
0.1687074064747948E+01 
0.3896282168946610E+01 

0.5378159157423945E-02 

0.1689195281391659E-01 

0.4443355769152601E-01 

0.1123335780703449E+00 

0.2711448390160103E+00 

0.6264263367286485E+00 

0.1408860943521751E+01 

0.3344043060866228E+01 

0.726E-04 

14 

0.1075073588251350E-02 

0.5889243490962496E-02 

0.1560078432135377E-01 

0.3258052212086110E-01 

0.6154351752779967E-01 

0.1109619891032348E+00 

0.1951651530857407E+00 

0.3377699882687942E+00 

0.577280541921 1481E+00 

0.9761165652290038E+00 

0.1635615445691163E+01 

0.2723809484786727E+01 

0.4541163041303490E+01 

0.7767616655342678E+01 

0.2783455121689438E-02 

0.7006395914900820E-02 

0.1279502133157069E-01 

0.2192733340131016E-01 

0.3737740049082059E-01 

0.6379243969367225E-01 

0 . 1084594588227473E+00 
0.1830223278438481E+00 
0.3061647832783700E+00 
0.5079755103629931E+00 
0.8381 174751258640E+00 
0.1385562498413431E+01 
0.2347348786059432E+01 
0.4444622409829190E+01 

0.366E-07 
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Table  2:  Quadratures  for  exponentials  (continued) 


Quadratures  for  the  integral 

fOO 

I  e~^^dx, 

Jo 

under  the  condition  that  1  <  t  <  500. 


m 

Nodes  (xj-) 

Weights  (in,) 

Error 

23 

0.6351980115825126E-03 

0.3390398468501349E-02 

0.8533988111011606E-02 

0.1642520313864894E-01 

0.2767253774178896E-01 

0.4324590958238922E-01 

0.6462964536146855E-01 

0.9401198766143719E-01 

0.1345150809637970E+00 

0.1904922248317389E+00 

0.2679364541057184E+00 

0.3750493014038291E+00 

0.5230239495194809E+00 

0.7271139907371750E+00 

0.1008092183746637E+01 

0.1394262856335610E+01 

0.1924307475579603E+01 

0.2651497602918851E+01 

0.3650430825 998876E+01 

0.5029133182526411E+01 

0.6954672288346456E+01 

0.9721470480335499E+01 

0.1402158019660932E+02 

0.1635014749191032E-02 
0.39060051 73541682E-02 
0.6439810761063304E-02 
0.9442176581726002E-02 
0.1321157735118120E-0a 
0.1817939708425140E-01 
0.2494901611240477E-01 
0.34331 18358535485E-01 
0.4739385869478771E-01 
0.6554991255546435E-01 
0.9069257914926255E-01 

0 . 12538 1 0727774845E+00 
0.1730865464562759E+00 
0.2385227411097447E+00 
0.3281047583403605E+00 
0.4506156991375293E+00 
0.6182455916131797E+00 
0.8483813375121229E+00 
0.1167112116649273E+01 
0.1617207388482339E+01 
0.2279901680035951E+01 
0.3352383978313540E+01 
0.5608355831510393E+01 

0.356E-12 
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Table  3:  Quadratures  for  exponentials  (continued) 


Quadratures  for  the  integral 

AOO 

/  e-^Ux, 

Jo 

under  the  condition  that  1  <  t  <  500. 


1  ^ 

Nodes  (x,) 

Weights  (wi) 

Error 

27 

0.5378759010624780E-03 

0.2860176825815242E-02 

0.714865861771 6300E-02 

0.1360965515937845E-01 

0.2257800188133212E-01 

0.3456421989535069E-01 

0.5032042618508775E-01 

0.7092509447124836E-01 

0.9788439120828463E-01 

0.1332509921950535E+00 

0.1797695570864978E+00 

0.2410654714132133E+00 

0.3218961915636380E+00 

0.4284852078938826E+00 

0.5689615509235298E+00 

0.7539347736933301E+00 

0.9972472224438443E+00 

0.1316964566299846E+01 

0.1736698582009859E+01 

0.2287418444638146E+01 

0.3010034073439038E+01 

0.3959315495048493E+01 

0.5210381702393131E+01 

0.6870768194824406E+01 

0.9106577764323245E+01 

0.1221294512896673E+02 

0.1689348652665484E+02 

0.1383311204046008E-02 

0.3279869733166365E-02 

0.5330932895600203E-02 

0.76460931 10803760E-02 

0.1037458793227033E-01 

0.1372178039022047E-01 

0.1796868836009351E-01 

0.2348971809947674E-01 

0.3076860552710760E-01 

0.4041894092839717E-01 

0.5321827718681367E-01 

0.7016094768858448E-01 

0.9253048536912244E-01 

0.1219928996130354E+00 

0.1607156476580828E+00 

0.2115215602167892E+00 

0.2780925850550500E+00 

0.3652478333806065E+00 

0.4793398853949993E+00 

0.6288554258416082E+00 

0.8254021 100491956E+00 

0.1085495633209734E+01 

0.1434174907278760E+01 

0.1913323186889750E+01 

0.2604342790201154E+01 

0.3708436699287805E+01 

0.60230861 566 1 5004E+0 1 

0.323E-14 
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Table  4:  Quadratures  for  complex  exponentials 


Quadratures  for  the  integral 

/■oo 

/  e-^^dx, 

Jo 

under  the  condition  that  z  €  C  lies  in  the  region  D  of  the  complex  plane  depicted  in  Figure  1. 


N 

Nodes  (xj) 

Weights  (wi) 

Error 

7 

0.1099271618238942E+00 

0.5491694162336780E+00 

0.1271416827286341E+01 

0.2239523056474301E+01 

0.3446836330005198E+01 

0.4877666068772302E+01 

0.6502607915187052E+01 

0.2775596224308371E+00 

0.5900612562744907E+00 

0.8478610527159362E+00 

0.1088510946164213E+01 

0.1323732065869006E+01 

0.1534838877513932E+01 

0.1719048349027934E+01 

0.107E-02 

10 

0.7940097370047949E-01 

0.4059967502704461E+00 

0.9586054827056690E+00 

0.1707633862341116E+01 

0.2634252243120157E+01 

0.3733067811454947E+01 

0.5005663556309191E+01 

0.6447614701968830E+01 

0.8049956086568744E+01 

0.9806270415536372E+01 

0.2021326824744206E+00 

0.4452920131070853E+00 

0.6549257007902238E+00 

0.8399190894283777E+60 

0.1012522786957398E+01 

0.1185698158021533E+01 

0.1358749093234873E+01 

0.1523775992304074E+01 

0.1681530325372958E+01 

0.1839363349513445E+01 

0.398E-04 

17 

0.4810701202067075E-01 

0.2505848757761927E+00 

0.6047137247359728E+00 

0. 1 0979569045692 1 7E+0 1 

0.1718338718562377E+01 

0.24561 16758149593E+01 

0.3304098340771468E+01 

0.4257638182548677E+01 

0.5314792420007071E+01 

0.6476281670685671E+01 

0.7744192070244406E+01 

0.9119761690912472E+01 

0.1060243019989989E+02 

0.1219107734541850E+02 

0.1388654365182865E+02 

0.1569478151303984E+02 

0.1 762993064234310E+02 

0.1231104634892695E+00 
0.2802783031579153E+00 
0.4258197681 747696E+00 
0.5586462093341786E+00 
0.6804909782597703E+00 
0.7938674951201829E+00 
0.9013184216898933E+00 
0.1005434048854524E+01 
0.1109016499341237E+01 
0.1214353975155841E+01 
0.1321691619491553E+01 
0.1429331461933850E+01 
0.1535781285935049E+01 
0.1641633199457229E+01 
0.1750408322216728E+01 
0.1869608472411435E+01 
0.2013038988665553E+01 

0.156E-07 

31 


Table  5:  Quadratures  for  complex  exponentials  (continued) 

Quadratures  for  the  integral 

/•oo 

/  e-^^dx, 

Jo 

under  the  condition  that  2  €  C  lies  in  the  region  D  of  the  complex  plane  depicted  in  Figure  1. 


N 

Nodes  (x,) 

Weights  (wi) 

Error 

26 

0.3186852812707167E-01 

0.1670696511352561E+00 

0.4070892206145096E+00 

0.74724391 75044331E+00 

0.1182020706329825E+01 

0.1705877918170689E+01 

0.2313727921958910E+01 

0.3001 165447553827E+01 

0.3764536084546370E+01 

0.4600933023898150E+01 

0.5508181675990506E+01 

0.6484841722444765E+01 

0.7530228834188133E+01 

0.8644426324529564E+01 

0.9828211695813383E+01 

0. 1 108280524626372E+02 

0.1240944939152817E+02 

0.1380901099033861E+02 

0.1528185219684248E+02 

0.1682806304875190E+02 

0.1844797932809189E+02 

0.2014289281830454E+02 

0.2191599603042421E+02 

0.2377386163950739E+02 

0.2572922307071097E+02 

0.2780603356526977E+02 

0.8168641985324435E-01 

0.1882692109574778E+00 

0.2909692905166523E+00 

0.3884080535348131E+00 

0.4802140096002388E+00 

0.5666489183411501E+00 

0.6483159527750231E+00 

0.7259502257876793E+00 

0.8003069049756682E+00 

0.8721245141158127E+00 

0.9421323128655350E+00 

0.1011073349889977E+01 

0.1079715704282156E+01 

0.1148813015862887E+01 

0.1218968645024915E+01 

0 . 1290432847974804E+01 
0.1363002295760370E+01 
0.1436177127644119E+01 
0.1509511398129948E+01 
0. 1582954591205827E+01 
0.1657078657112459E+01 
0.1733260266399871E+01 
0.1814015981921987E+01 
0.1903865269011978E+01 
0.201 1585306477846E+01 
0.2156295623247961E+01 

0.801E-12 

i 

32 


Table  6:  Quadratures  for  complex  exponentials  (continued) 


Quadratures  for  the  integral 

too 

/  e-^^dx, 

Jo 

under  the  condition  that  z  £  C  lies  in  the  region  D  of  the  complex  plane  depicted  in  Figure  1. 


N 

Nodes  (xi) 

Weights 

Error 

32 

0.2599836293936463E-01 

0.1365267137471029E+00 

0.3335929178284293E+00 

0.6144848200204676E+00 

0.9757788966059570E+00 

0.1413764904033582E+01 

0.1924760915837073E+01 

0.2505322024543478E+01 

0.3152355883303843E+01 

0.3863171917607857E+01 

0.4635491817809779E+01 

0.5467443320214809E+01 

0.6357552049608270E+01 

0.7304739365689955E+01 

0.8308328039815608E+01 

0.9368051427922781E+01 

0.1048405410784282E+02 

0.1165686219043638E+02 

0.1288729607812993E+02 

0.1417631366171486E+02 

0.1552481519724356E+02 

0 . 1 693348466077443E+02 
0.1840274326790348E+02 
0.1993284537567664E+02 
0.2152409758022736E+02 
0.2317716791120881E+02 
0.2489347969157359E+02 
0.2667574373640977E+02 
0.2852878296929011E+02 
0.3046102218822619E+02 
0.3248756039307094E+02 
0.3463653566705793E+02 

0.6666723984893712E-01 

0.1541635906448656E+00 

0.2395065808061905E+00 

0.3216997252047818E+00 

0.4002647717267098E+00 

0.4750901483689786E+00 

0.5463258676163205E+00 

0.6142802064945914E+00 

0.6793389599054632E+00 

0.7419135289744640E+00 

0.8024143819948454E+00 

0.86124293381 13290E+00 

0.9187946472618356E+00 

0.9754669332016366E+00 

0.1031665911277885E+01 

0.1087805382732233E+01 

0.1144289073127857E+01 

0.1201466511491073E+01 

0.1259564363290567E+01 

0.1318623693719137E+01 

0.1378494220804057E+01 

0.1438913247537656E+01 

0.1499642611759877E+01 

0.1560607527964651E+01 

0.1622001781246697E+01 

0.1684367617089771E+01 

0.1748695093705612E+01 

0.1816615869561691E+01 

0.1890833030920412E+01 

0.1976151351204436E+01 

0.2082162437141457E+01 

0.2231151873780132E+01 

0.282E-14 

33 


Table  7:  Quadratures  for  exponentials  midtiplied  by  Iq 


Quadratures  for  tbe  integral 

f  Io{xy)e~^^dx, 

Jo 

under  the  condition  that  t  €  [1,500]  and  y  6  [0,t  -  1]. 


N 

Nodes  (xt) 

Weights  {w{) 

Error 

6 

0.2014951814758335E-02 

0.1524995123116156E-01 

0.6752356088881507E-01 

0.2501565900307728E+00 

0.8236005667135590E+00 

0.2511917343393000E+01 

0.5609959210156781E-02 

0.2486640535654579E-01 

0.9333052982723889E-01 

0.3114318476382451E+00 

0.9419185523143980E+00 

0.2783935122763889E+01 

0.997E-03 

8 

0.1310754453518395E-02 

0.8427815046421337E-02 

0.2934337237542595E-01 

0.8663331921085451E-01 

0.2362895983874226E+00 

0.6083828926140039E+00 

0.1496074584403707E+01 

0.3613104405570935E+01 

0.3522290544269296E-02 

0.1188533050797089E-01 

0.3342585308732757E-01 

0.8985909963205312E-01 

0.2297963615059060E+00 

0.5600968338362269E+00 

0.1320244303895297E+01 

0.3253613362413727E+01 

0.892E-04 

14 

0.6424288534795956E-03 

0.3562319666990144E-02 

0.9643424057074440E-02 

0.2074298599770349E-01 

0.4057928260022333E-01 

0.7600572280169251E-01 

0.1390443053485344E+00 

0.2503136051566992E+00 

0.4447622918282108E+00 

0.7811276346003586E+00 

0.1357818162257100E+01 

0.2341992534236977E+01 

0.4036529413075654E+01 

0.7126502974662635E+01 

0.1667964367860395E-02 

0.4298903067080389E-02 

0.8159545461265918E-02 

0.1463864640961027E-01 

0.2614391453322226E-01 

0.4665755537868725E-01 

0.8276361628521883E-01 

0.1454478222995341E+00 

0.2529871458046016E+00 

0.4357009925372973E+00 

0 . 7446059596729966E+00 
0.1271434906786924E+01 
0.2216807353831690E+01 
0.4302367103374836E+01 

0.900E-07 
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Table  8:  Quadratures  for  exponentials  multiplied  by  Iq  (continued) 


Quadratures  for  the  integral 

roo 

/  Io{xy)e~'^^dx, 

Jo 

under  the  condition  that  t  €  [1,500]  and  j/  €  [0,t  -  Ij. 


N 

Nodes  (a:,) 

Weights  (wi) 

Error 

24 

0.3495789092315762E-03 

0.1870021782900649E-02 

0.4726267824851739E-02 

0.9151716313193496E-02 

0.1554236980028004E-01 

0.2452696840462340E-01 

0.3706303150690139E-01 

0.5456277459785074E-01 

0.7905495952668643E-01 

0.1134006713086074E+00 

0. 1 615933468425095E+00 

0.2291 729794308585E+00 

0.3238013639554152E+00 

0.4560558428504051E+00 

0.64051 13416979851E+00 

0.8972279324556263E+00 

0.1253809582199313E+01 

0.1748318313157252E+01 

0.2433578887375728E+01 

0.3384010828873551E+01 

0.4707795467897619E+01 

0.6572440958883276E+01 

0.9272736577094724E+01 

0.1349943616527142E+02 

0.90028406k465873E-03 

0.2160252603038041E-02 

0.3590338301868660E-02 

0.5325399369607390E-02 

0.7559106650189445E-02 

0.1056758825931061E-0i 

0.1473724813105991E-01 

0.2059685692839059E-01 

0.2885941258668471E-01 

0.4049079881256631E-01 

0.5680844261155735E-01 

0.7961469039423471E-01 

0.1113872397913262E+00 

0.1555264308649468E+00 

0.2166942382013035E+00 

0.3012922008808448E+00 

0.4181453830530776E+00 

0.5795680006227947E+00 

0.8031688787093632E+00 

0.1115348337169697E+01 

0.1559213913124744E+01 

0.2216187240183199E+01 

0.3282974439738362E+01 

0.5529398603135539E+01 

0.925E-12 
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Table  9:  Quadratures  for  exponentials  multiplied  by  /q  (continued) 


Quadratures  for  the  integral 

too 

/  /o(xy)e“®‘dx, 

Jo 

under  the  condition  that  t  €  [1,500]  and  y  £  [0,t  —  Ij. 


N 

Nodes  (xi) 

Weights  (wi) 

Error 

29 

0.2855179413353365E-03 

0.1519624696728258E-02 

0.3804359141657344E-02 

0.7260138000706486E-02 

0.1208205371062810E-01 

0.1856564543199398E-01 

0.2714156753309568E-01 

0.384201 7800878239E-01 

0.5324783256625659E-01 

0.7277755829761968E-01 

0.9855788611173273E-01 

0.1326465035778468E+00 

0.1777590387840778E+00 

0.2374657658898870E+00 

0.3164509240422835E+00 

0.4208524457939620E+00 

0.5587051648321601E+00 

0 .7405 1 85479404663E+00 
0.9800319873390735E+00 
0.1295209795621391E+01 
0.1709570851677607E+01 
0.2254009385987865E+01 
0.2969389638669206E+01 
0.3910476327629495E+01 
0.5152430007642100E+01 
0.6802867813529709E+01 
0.9027979519502084E+01 
0.1212289908066820E+02 
0.1679085599535762E+02 

0.7344503079351386E-03 

0.1744538390662211E-02 

0.2844687196642974E-02 

0.4098961298933580E-02 

0.5593550200298448E-02 

0.7444670271885530E-02 

0.9807968524698940E-02 

0.1288914176031762E-01 

0.1695687345717790E-01 

0.2235879759838917E-01 

0.2954235698585380E-01 

0.3908423859367898E-01 

0.5173159700695577E-01 

0.6845695550893067E-01 

0.9052903520482935E-01 

0.1196036182345700E+00 

0.1578409524693449E+00 

0.2080593451794129E+00 

0.2739397750144418E+00 

0.3603059290242059E+00 

0.4735231867763476E+00 

0.6221016600956893E+00 

0.8176841100086656E+00 

0.1076831175000068E+01 

0.1424628439002124E+01 

0.1902988149814232E+01 

0.2593285548365225E+01 

0.3696550722303479E+01 

0.6009492062220468E+01 

0.299E-14 
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Table  10:  Quadratures  for  exponentials  multiplied  by  Jq 


Quadratures  for  tbe  integral 

roo 

I  Jo{xy)e~^*dx, 

Jo 

under  the  condition  that  t  £  [1,4]  and  y  €  [0,4\^]. 


N 

Nodes  (ij) 

Weights  (in,) 

Error 

8 

0.1093474676900044E+00 

0.5176974101534121E+00 

0.1133065916111916E+01 

0.1881350151107404E+01 

0.2717854096012053E+01 

0.3616502749074490E+01 

0.4562710533038212E+01 

0.5549008853485283E+01 

0.2710750266277354E+00 

0.5276915884394641E+00 

0.6915150441387948E+00 

0.7983440040645204E+00 

0.8716416012135397E+00 

0.9264383911692414E+00 

0.9729462225948307E+00 

0.1024138658446855E+01 

0.162E-02 

12 

0.7685522448236467E-01 

0.3802271685596512E+00 

0.8629501667245919E+00 

0.1477406574242533E+01 

0.2190593072512602E+01 

0.2979188555054684E+01 

0.3826805213168235E+01 

0.4722181214285143E+01 

0.5657828852278510E+01 

0.6629008403962641E+01 

0.763291 1519449263E+01 

0.8669258567695921E+01 

0.1937803229242497E+00 

0.4024780894501363E+00 

0.5551232854865536E+00 

0.6684012296815303E+00 

0.7541446224405415E+00 

0.8203905361353097E+00 

0.8731017778158731E+00 

0.9167109597437153E+00 

0.9545728875259875E+00 

0.9893709749159459E+00 

0.1023874368056413E+01 

0.1067824933823433E+01 

0.709E-04 

37 
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Table  11:  Quadratures  for  exponentials  multiplied  by  Jo  (continued) 


Quadratures  for  the  integral 

roo 

/  Jo{xy)e~'^*dx, 

Jo 

under  the  condition  that  t  €  [1,4]  and  y  €  [0,4v/5]. 


N 

Nodes  (ij) 

Weights  (wi) 

Error 

21 

0.45571 10658309451E-01 

0.2345692419777160E+00 

0.5560507435597863E+00 

0.9888622621 190326E+00 

0.1514051750681985E+01 

0.2116308405894669E+01 

0.2783466404955423E+01 

0.3505774922334249E+0 1 

0.4275358227777114E+01 

0.5085850421343891E+01 

0.5932129958898720E+01 

0.6810110232637652E+01 

0.7716569940856932E+01 

0.8649018954485772E+01 

0.9605599035641322E+01 

0.1058501999437337E+02 

0.1 158653269800637E+02 

0.1260993904708917E+02 

0.1365563090220990E+02 

0.1472471197434301E+02 

0.1582111587898742E+02 

0.1162693279863745E+00 
0.2587357630437822E+00 
0.3805968474435821E+00 
0.4818526125575090E+00 
0.5659864564983776E+00 
0.6365120448492290E+00 
0.6961876640755223E+00 
0.7471226714069135E+00 
0. 7909873646366727E+00 
0.8291470860681295E+00 
0.8627341849390908E+00 
0.8926928974293094E+00 
0.9198161427100997E+00 
0.944780931241 1800E+00 
0.9681847744461964E+00 
0.9905857720543664E+00 
0.1012551797552972E+01 
0.1034729712476950E+01 
0.1057973242328095E+01 
0.1083960111883219E+01 
0.1123223099240344E+01 

0.553E-07 

38 


Table  12:  Quadratures  for  exponentials  multiplied  by  Jo  (continued) 

Quadratures  for  the  integral 

f  Jo{xy)e~^^dx, 

Jo 

under  the  condition  that  t  €  [1,4]  and  y  6  [0,4y/2]. 


N 

Nodes  (xi) 

Weights  (wi) 

Error 

31 

0.3135427831034307E-01 

0.8024339887513055E-01 

0.195E-10 

0.1633283571233953E+00 

0.1826480678501762E+00 

0.3938519939248281E+00 

0.2767843166469057E+00 

0.7134977521472219E+00 

0.3607878783606646E+00 

0.1112086865038666E+01 

0.4347972857765461E+00 

0.1580107134125432E+01 

0.4998566920961263E+00 

0.2109237828913374E+01 

0.5572245708288816E+00 

0.2692380266283717E+01 

0.6080577723352450E+00 

0.3323499813712884E+01 

0.6533243161908575E+00 

0.3997438592537748E+01 

0.6938153813251570E+00 

0.4709757745640057E+01 

0.7301846396770468E+00 

0.5456619026749223E+01 

0 . 7629849054856576E+00 

» 

0.6234697842382201E+01 

0.7926941021324604E+00 

0.704 1 1 1 7206272996E+0 1 

0.8197316250932437E+00 

0.7873394668593335E+01 

0.8444685010430712E+00 

0.8729397582414611E+01 

0.8672343080184358E+00 

0.9607304495127067E+01 

0.8883227316567794E+00 

0.1050557184595927E+02 

0.9079967283507420E+00 

0.1142290582410390E+02 

0.9264937318864332E+00 

0.1235823954272516E+02 

0.9440311133071873E+00 

0.1331071586594676E+02 

0.9608120820932964E+00 

0.1427967645292451E+02 

0.9770323352089769E+00 

0.1526465797627137E+02 

0.9928880149352948E+00 

0.1626539716056547E+02 

0.1008585996949826E+01 

0.1728184746281412E+02 

0.1024358393992377E+01 

0.1831421215592517E+02 

0.1040484810479409E+01 

0.1936300165981823E+02 

0.1057329638354610E+01 

0.2042912676709522E+02 

0.1075408955912917E+01 

0.2151404706129260E+02 

0. 1095542344669040E+01 

0.2262012600158427E+02 

0.1119613828562103E+01 

0.2375360482790972E+02 

0.1159632125663126E+01 

39 
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Table  13:  Quadratures  for  exponentials  multiplied  by  Jq  (continued) 


Quadratures  for  the  integral 

/•oo 

/  Jo{xy)e~’‘^dx, 

Jo 

under  the  condition  that  t  €  [1,4]  and  y  €  [0,4\/2]. 


N 

Nodes  {xi) 

Weights  (wi) 

Error 

40 

0.2450466782681923E-01 
0.1282174857176023E+00 
0.3113726824539808E+00 
0.5689058522320320E+00 
0.8947716487384897E+00 
0.1282740805986842E+01 
0.1726889168240672E+01 
0.2221829818237509E+01 
0.2762775128318641E+01 
0.3345508771449190E+01 
0.3966321542990838E+01 
0.4621940741778159E+01 
0.5309466588865855E+01 
0.6026319872772099E+01 
0.6770200467500875E+01 
0.7539054647612656E+01 
0.8331048881923509E+01 
0.9144548183658569E+01 
0.9978097658074295E+01 
0.1083040639472408E+02 
0.1170033322583332E+02 
0.1258687411609641E+02 
0.1348915109370842E+02 
0.1440640271349045E+02 
0.1533797609080401E+02 
(Continued  in 

0.6278240289055225E-01 

0.1441452770046284E+00 

0.2213018947962683E+00 

0.2927436383298983E+00 

0 .357941 71 72934986E+00 
0.4170055991300023E+00 
0.4703924681663025E+00 
0.5186914865873826E+00 
0.5624975493703959E+00 
0.6023527423690355E+00 
0.6387284183350381E+00 
0.6720275534119337E+00 
0.7025949580565075E+00 
0.7307287944249567E+00 
0.7566905396059194E+00 
0.7807125897154367E+00 
0.8030036866777486E+00 
0.8237527017133555E+00 
0.8431313322224179E+00 
0.8612961564381806E+00 
0.8783903510592598E+00 
0.894545259933 1251E+00 
0.9098819208500799E+00 
0.9245126105589752E+00 
0.9385424474160881E+00 
the  next  table) 
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Table  14:  Quadratures  for  exponentials  multiplied  by  Jq  (continued) 


Quadratures  for  tbe  integral 

f  Jo{xy)e~^^dx, 

Jo 


under  the  condition  that  t  6  [1,4]  and  y  6  [0,4-\/2]. 


N 

Nodes  (xi)  Weights 

Error 

40 

(Continued  from  t! 
0.1628332058271154E+02 
0.1724198323818701E+02 
0.1821360621166772E4-02 
0.1919792646634440E+02 
0.2019477834484074E+02 
0.2120410011858627E+02 
0.2222594665494397E+02 
0.2326051 184537432E+02 
0.2430816659609159E+02 
0.2536951 778065644E+02 
0.2644549219792192E+02 
0.2753746108815434E+02 
0.2864745668803461E+02 
0.2977833401798621E+02 
0.3093837103779525E+02 

le  preceding  table) 
0.9520710907702972E+00 
0.9651945915710470E+00 
0.9780074909877168E+00 
0.9906053454228810E+00 

0 . 1 0030880346 1 0363E+01 
0.1015564634343247E+01 
0.1028161001682954E+01 
0.1041032846271463E+01 
0.1054383431852286E+01 
0.1068492599609564E+01 
0.1083742894323642E+01 
0.1100713259554107E+01 
0.1120408914213389E+01 
0.1144615248532595E+01 
0.1182108938213342E+01 

0.147E-13 
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